{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "## Quantitative phase analysis (QPA)\n", "Fox/ObjCryst++ was not designed with QPA in mind, but it\n", "is still possible to do it when the phases are known\n", "and the profiles not too complicated.\n", "\n", "Here we just try the 'simple' cases of the QPA Round Robin of\n", "1999 (https://www.iucr.org/__data/iucr/powder/QARR/index.html)" ] }, { "cell_type": "code", "execution_count": 1, "id": "1", "metadata": {}, "outputs": [], "source": [ "# 'widget' allows live update and works in both classical notebooks and jupyter-lab.\n", "# Otherwise 'notebook', 'ipympl', 'inline' can be used\n", "%matplotlib widget\n", "\n", "import os\n", "import pyobjcryst\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pyobjcryst.crystal import *\n", "from pyobjcryst.powderpattern import *\n", "from pyobjcryst.indexing import *\n", "from pyobjcryst.molecule import *\n", "from pyobjcryst.globaloptim import MonteCarlo\n", "from pyobjcryst.io import xml_cryst_file_save_global\n", "from pyobjcryst.lsq import LSQ\n", "from pyobjcryst.refinableobj import refpartype_scattdata_scale\n" ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "### Data and CIF sources" ] }, { "cell_type": "code", "execution_count": 2, "id": "3", "metadata": {}, "outputs": [], "source": [ "# Data from QPA round-robin (cpd-1a.prn):\n", "# https://www.iucr.org/__data/iucr/powder/QARR/col/cpd-1a.prn\n", "# Try samples 1a to 1h\n", "\n", "data_file = \"data/cpd-1a.prn\"\n", "\n", "# Data from Crystal structures from the Crystallography Open Database:\n", "# https://www.crystallography.net/cod/\n", "\n", "cod_phases = [1000032, 9008877, 5000222] # Al2O3, ZnO, CaF2 - needs checking\n" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "### Create Powder pattern" ] }, { "cell_type": "code", "execution_count": 3, "id": "5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.5418366591135662\n", "Imported powder pattern: 7251 points, 2theta= 5.000 -> 150.000, step= 0.020\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "e982d5d79ef84be0a9192ebafce0ea53", "version_major": 2, "version_minor": 0 }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAGQCAYAAAD2lq6fAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAATgVJREFUeJzt3Xl4VPXd9/HPZF9IJhvJMBAgFGQLVQSLIBUqsqhIbX2qFhtta3EHoyhq1Uq9CrhU9FZaVB5vtSLFu33EuiKIFkTWBlFAVg0QhBCEZLJPtvP8gXPuLCwxTDIz57xf15XrImd+mfM7X7LMZ77n/I7DMAxDAAAAAADbCQv0BAAAAAAAgUEgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsKmIQE8AZ6ahoUEHDx5UQkKCHA5HoKcDAAAAnDHDMFRWVia3262wMHpY7YlAGOIOHjyozMzMQE8DAAAA8LuCggJ169Yt0NOwNAJhiEtISJB0/IclMTExwLMBAAAAzlxpaakyMzPN17poPwTCEOc7TTQxMZFACAAAAEvhkqj2xwm5AAAAAGBTBEIAAAAAsCkCIQAAAADYFNcQAgAAAH5SX1+v2traQE8j6EVGRio8PDzQ04AIhAAAAMAZMwxDhYWFKikpCfRUQkZSUpJcLhcLxwQYgRAAAAA4Q74wmJ6erri4OELOKRiGocrKShUVFUmSunTpEuAZ2RuBEAAAADgD9fX1ZhhMTU0N9HRCQmxsrCSpqKhI6enpnD4aQCwqAwAAAJwB3zWDcXFxAZ5JaPHVi2suA4tACAAAAPgBp4l+P9QrOBAIAQAAAMCmCISAxfzrX//SG2+8EehpAACAEPfvf/9bDoeDlVMtjkVlAIu54oorJB1fwQsAAAA4FTqEAAAAAGBTBEIAAADAprxer6ZNm6b09HTFxMRo5MiR2rhxY5Mxn376qc4++2zFxMRo2LBh2rJli/nYvn37dPnllys5OVnx8fEaOHCg3nvvvY4+DJwB2wXCVatW6fLLL5fb7ZbD4dCbb77Z5HHDMDRz5ky53W7FxsZq9OjR2rZtW5MxXq9XU6dOVVpamuLj4zVp0iQdOHCgyZji4mLl5OTI6XTK6XQqJyenxfnX+/fv1+WXX674+HilpaVp2rRpqqmpaY/DBgAAAFqYMWOG/t//+3965ZVXtGnTJvXu3Vvjx4/XsWPHzDH33HOP/vznP2vjxo1KT0/XpEmTzFtF3HbbbfJ6vVq1apW2bNmixx57TJ06dQrU4aANbHcNYUVFhc4++2z95je/0ZVXXtni8ccff1xz587Vyy+/rLPOOkt/+tOfNHbsWO3cuVMJCQmSpNzcXL399ttavHixUlNTNX36dE2cOFF5eXnmTTUnT56sAwcOaOnSpZKkG2+8UTk5OXr77bclHb+B6WWXXabOnTtr9erVOnr0qK6//noZhqFnn322g6oBAACA9lJZWakdO3Z0+H779evXqnsiVlRUaP78+Xr55Zd1ySWXSJIWLFig5cuX68UXX9R5550nSXr44Yc1duxYSdIrr7yibt26acmSJbrqqqu0f/9+XXnllRo0aJAkqVevXu10VGg3ho1JMpYsWWJ+3tDQYLhcLuPRRx81t1VXVxtOp9N47rnnDMMwjJKSEiMyMtJYvHixOeabb74xwsLCjKVLlxqGYRhffvmlIclYt26dOWbt2rWGJGPHjh2GYRjGe++9Z4SFhRnffPONOebvf/+7ER0dbXg8nlYfg8fjMSR9r6+BtUkybP6j7Re1tbVGfX19oKcBAAgBVVVVxpdffmlUVVU12Z6Xl2f+Xe7Ij7y8vFbN+/PPPzckGXv37m2y/YorrjB+85vfGB9//LEhydi3b1+Tx8855xxj5syZhmEYxoIFC4yIiAhjxIgRxh/+8Afj888/P+O6GQavcTuS7TqEp5Kfn6/CwkKNGzfO3BYdHa1Ro0ZpzZo1uummm5SXl6fa2tomY9xut7Kzs7VmzRqNHz9ea9euldPp1LBhw8wx559/vpxOp9asWaO+fftq7dq1ys7OltvtNseMHz9eXq9XeXl5+slPfnLCOXq9Xnm9XvPz0tJSf5YAwHciIyN1ySWXcB0EAKDN+vXrp7y8vIDstzWM71Ykb36DeMMwTnvTeN/jv/vd7zR+/Hi9++67WrZsmebMmaMnn3xSU6dObcPMEQgEwkYKCwslSRkZGU22Z2RkaN++feaYqKgoJScntxjj+/rCwkKlp6e3eP709PQmY5rvJzk5WVFRUeaYE5kzZ47++Mc/fs8jA9AW77//fqCnAAAIYXFxcTr33HMDPY2T6t27t6KiorR69WpNnjxZklRbW6v//Oc/ys3NNcetW7dO3bt3l3R8nYxdu3Y1CZ2ZmZm6+eabdfPNN+v+++/XggULCIQhxHaLyrRGW94laT7mROPbMqa5+++/Xx6Px/woKCg45bwAAACAE4mPj9ctt9yie+65R0uXLtWXX36pKVOmqLKyUjfccIM57pFHHtGKFSu0detW/frXv1ZaWpp53+Pc3Fx98MEHys/P16ZNm/TRRx+pf//+AToitAUdwkZcLpek4927Ll26mNuLiorMbp7L5VJNTY2Ki4ubdAmLioo0YsQIc8zhw4dbPP+RI0eaPM/69eubPF5cXKza2toWncPGoqOjFR0d3cYjBAAAAP7Xo48+qoaGBuXk5KisrExDhw7VBx980OR17qOPPqo77rhDu3fv1tlnn6233npLUVFRko4vlHjbbbfpwIEDSkxM1IQJE/TUU08F6nDQBnQIG8nKypLL5dLy5cvNbTU1NVq5cqUZ9oYMGaLIyMgmYw4dOqStW7eaY4YPHy6Px6MNGzaYY9avXy+Px9NkzNatW3Xo0CFzzLJlyxQdHa0hQ4a063ECAAAAkhQTE6NnnnlGR44cUXV1tVavXm2uLjp69GgZhqGJEydq69at8nq92rBhg84++2zz65999lnt2bNH1dXVKioq0t/+9jelpqYG6nDQBrbrEJaXl2vPnj3m5/n5+dq8ebNSUlLUvXt35ebmavbs2erTp4/69Omj2bNnKy4uzjyv2ul06oYbbtD06dOVmpqqlJQU3X333Ro0aJAuvvhiSVL//v01YcIETZkyRc8//7yk47edmDhxovr27StJGjdunAYMGKCcnBw98cQTOnbsmO6++25NmTJFiYmJHVwVAAAAAHZku0D4n//8p8kKnnfddZck6frrr9fLL7+sGTNmqKqqSrfeequKi4s1bNgwLVu2zLwHoSQ99dRTioiI0FVXXaWqqiqNGTNGL7/8snkPQkl67bXXNG3aNHM10kmTJmnevHnm4+Hh4Xr33Xd166236oILLlBsbKwmT56sP//5z+1dAgAAAACQJDkM33qzCEmlpaVyOp3yeDx0FiHpfxcr4kf7zFBHAEBrVVdXKz8/X1lZWYqJiQn0dELGqerGa9yOwzWEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAPC9zZw5U+ecc06gp4EzRCAEAAAAAJsiEAIAAACATREIAQAAAJtqaGjQY489pt69eys6Olrdu3fXrFmzJEn33nuvzjrrLMXFxalXr1566KGHVFtbe8rn++///m8NHDhQ0dHR6tKli26//faOOAycgYhATwAAAABAYNx///1asGCBnnrqKY0cOVKHDh3Sjh07JEkJCQl6+eWX5Xa7tWXLFk2ZMkUJCQmaMWPGCZ9r/vz5uuuuu/Too4/qkksukcfj0aefftqRh4M2cBiGYQR6Emi70tJSOZ1OeTweJSYmBno6CAIOh0OSxI/2maGOAIDWqq6uVn5+vrKyshQTExPo6bRaWVmZOnfurHnz5ul3v/vdacc/8cQTev311/Wf//xH0vFFZd58801t3rxZktS1a1f95je/0Z/+9KdW7f9UdeM1bsehQwgAAAC0h8pK6btuW4fq10+KizvtsO3bt8vr9WrMmDEnfPyf//ynnn76ae3Zs0fl5eWqq6s7aTgrKirSwYMHT/pcCF4EQgAAAKA97NghDRnS8fvNy5POPfe0w2JjY0/62Lp163TNNdfoj3/8o8aPHy+n06nFixfrySef/N7PheBGIAQAAADaQ79+x8NZIPbbCn369FFsbKxWrFjR4pTRTz/9VD169NADDzxgbtu3b99JnyshIUE9e/bUihUr9JOf/KRt80ZAEAgBAACA9hAX16pOXaDExMTo3nvv1YwZMxQVFaULLrhAR44c0bZt29S7d2/t379fixcv1nnnnad3331XS5YsOeXzzZw5UzfffLPS09N1ySWXqKysTJ9++qmmTp3aQUeEtiAQAgAAADb10EMPKSIiQn/4wx908OBBdenSRTfffLNuuOEG3Xnnnbr99tvl9Xp12WWX6aGHHtLMmTNP+lzXX3+9qqur9dRTT+nuu+9WWlqa/s//+T8ddzBoE1YZDXGswITmWB3TP6gjAKC1QnWV0UBjldHgwI3pAQAAAMCmCIQAAAAAYFMEQgAAAACwKQIhAAAAANgUgRAAAAAAbIpACAAAAPgBK1N/P9QrOBAIAQAAgDMQGRkpSaqsrAzwTEKLr16++iEwuDE9AAAAcAbCw8OVlJSkoqIiSVJcXJx5P1u0ZBiGKisrVVRUpKSkJIWHhwd6SrZGIAQAAADOkMvlkiQzFOL0kpKSzLohcAiEAAAAwBlyOBzq0qWL0tPTVVtbG+jpBL3IyEg6g0GCQAgAAAD4SXh4OEEHIYVFZQAAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhA2U1dXpwcffFBZWVmKjY1Vr1699Mgjj6ihocEcYxiGZs6cKbfbrdjYWI0ePVrbtm1r8jxer1dTp05VWlqa4uPjNWnSJB04cKDJmOLiYuXk5MjpdMrpdConJ0clJSUdcZgAAAAAQCBs7rHHHtNzzz2nefPmafv27Xr88cf1xBNP6NlnnzXHPP7445o7d67mzZunjRs3yuVyaezYsSorKzPH5ObmasmSJVq8eLFWr16t8vJyTZw4UfX19eaYyZMna/PmzVq6dKmWLl2qzZs3Kycnp0OPFwAAAIB9OQzDMAI9iWAyceJEZWRk6MUXXzS3XXnllYqLi9Orr74qwzDkdruVm5ure++9V9LxbmBGRoYee+wx3XTTTfJ4POrcubNeffVVXX311ZKkgwcPKjMzU++9957Gjx+v7du3a8CAAVq3bp2GDRsmSVq3bp2GDx+uHTt2qG/fvq2ab2lpqZxOpzwejxITE/1cDYQih8Mh6XgnG21HHQEACBxe43YcOoTNjBw5UitWrNCuXbskSZ9//rlWr16tSy+9VJKUn5+vwsJCjRs3zvya6OhojRo1SmvWrJEk5eXlqba2tskYt9ut7Oxsc8zatWvldDrNMChJ559/vpxOpzkGAAAAANpTRKAnEGzuvfdeeTwe9evXT+Hh4aqvr9esWbP0y1/+UpJUWFgoScrIyGjydRkZGdq3b585JioqSsnJyS3G+L6+sLBQ6enpLfafnp5ujjkRr9crr9drfl5aWtqGowQAAAAAOoQtvP7661q4cKEWLVqkTZs26ZVXXtGf//xnvfLKK03G+U4n8zEMo8W25pqPOdH40z3PnDlzzEVonE6nMjMzW3NYAAAAANACgbCZe+65R/fdd5+uueYaDRo0SDk5Obrzzjs1Z84cSZLL5ZKkFl28oqIis2vocrlUU1Oj4uLiU445fPhwi/0fOXKkRfexsfvvv18ej8f8KCgoaPvBAgAAALA1AmEzlZWVCgtrWpbw8HDzthNZWVlyuVxavny5+XhNTY1WrlypESNGSJKGDBmiyMjIJmMOHTqkrVu3mmOGDx8uj8ejDRs2mGPWr18vj8djjjmR6OhoJSYmNvkAAAAAgLbgGsJmLr/8cs2aNUvdu3fXwIED9dlnn2nu3Ln67W9/K+n4aZ65ubmaPXu2+vTpoz59+mj27NmKi4vT5MmTJUlOp1M33HCDpk+frtTUVKWkpOjuu+/WoEGDdPHFF0uS+vfvrwkTJmjKlCl6/vnnJUk33nijJk6c2OoVRgEAAADgTBAIm3n22Wf10EMP6dZbb1VRUZHcbrduuukm/eEPfzDHzJgxQ1VVVbr11ltVXFysYcOGadmyZUpISDDHPPXUU4qIiNBVV12lqqoqjRkzRi+//LLCw8PNMa+99pqmTZtmrkY6adIkzZs3r+MOFgAAAICtcR/CEMc9WtAc98/zD+oIAEDg8Bq343ANIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIQAAAADYFIEQAAAAAGyKQAgAAAAANkUgBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIAAAAAsCkCIWBRhmEEegoAAAAIcgRCwKIIhAAAADgdAiFgUQRCAAAAnA6BEAAAAABsikAIWBQdQgAAAJwOgRCwKAIhAAAATodACAAAAAA2RSAELIoOIQAAAE6HQAhYFIEQAAAAp0MgBCyKQAgAAIDTIRACAAAAgE0RCAGLokMIAACA0yEQAhZFIGw7agcAAOyCQAgAAAAANkUgBCyKLlfbUTsAAGAXBELAogg1AAAAOB0CIQA0Q5gGAAB2QSAELIpQAwAAgNMhEAIWRSBsO2oHAADsgkAIWBShBgAAAKdDIDyBb775Rr/61a+UmpqquLg4nXPOOcrLyzMfNwxDM2fOlNvtVmxsrEaPHq1t27Y1eQ6v16upU6cqLS1N8fHxmjRpkg4cONBkTHFxsXJycuR0OuV0OpWTk6OSkpKOOEQAp0CYBgAAdkEgbKa4uFgXXHCBIiMj9f777+vLL7/Uk08+qaSkJHPM448/rrlz52revHnauHGjXC6Xxo4dq7KyMnNMbm6ulixZosWLF2v16tUqLy/XxIkTVV9fb46ZPHmyNm/erKVLl2rp0qXavHmzcnJyOvJwYWGEGgAAAJyOw+BVYxP33XefPv30U33yyScnfNwwDLndbuXm5uree++VdLwbmJGRoccee0w33XSTPB6POnfurFdffVVXX321JOngwYPKzMzUe++9p/Hjx2v79u0aMGCA1q1bp2HDhkmS1q1bp+HDh2vHjh3q27dvq+ZbWloqp9Mpj8ejxMREP1QAoc7hcEiSjh07puTk5ADPJjTV1NQoOjpaEsEaAIBA4DVux6FD2Mxbb72loUOH6he/+IXS09M1ePBgLViwwHw8Pz9fhYWFGjdunLktOjpao0aN0po1ayRJeXl5qq2tbTLG7XYrOzvbHLN27Vo5nU4zDErS+eefL6fTaY45Ea/Xq9LS0iYfAAAAANAWBMJmvv76a82fP199+vTRBx98oJtvvlnTpk3T3/72N0lSYWGhJCkjI6PJ12VkZJiPFRYWKioqqkV3pvmY9PT0FvtPT083x5zInDlzzGsOnU6nMjMz236wsDQ6W21H7QAAgF0QCJtpaGjQueeeq9mzZ2vw4MG66aabNGXKFM2fP7/JON9peT6GYbTY1lzzMScaf7rnuf/+++XxeMyPgoKC1hwWbIhQAwAAgNMhEDbTpUsXDRgwoMm2/v37a//+/ZIkl8slSS26eEVFRWbX0OVyqaamRsXFxaccc/jw4Rb7P3LkSIvuY2PR0dFKTExs8gGcCIEQAAAAp0MgbOaCCy7Qzp07m2zbtWuXevToIUnKysqSy+XS8uXLzcdramq0cuVKjRgxQpI0ZMgQRUZGNhlz6NAhbd261RwzfPhweTwebdiwwRyzfv16eTwecwyAwCBMAwAAu4gI9ASCzZ133qkRI0Zo9uzZuuqqq7Rhwwa98MILeuGFFyQdP80zNzdXs2fPVp8+fdSnTx/Nnj1bcXFxmjx5siTJ6XTqhhtu0PTp05WamqqUlBTdfffdGjRokC6++GJJx7uOEyZM0JQpU/T8889Lkm688UZNnDix1SuMAqdCqAEAAMDpEAibOe+887RkyRLdf//9euSRR5SVlaWnn35a1157rTlmxowZqqqq0q233qri4mINGzZMy5YtU0JCgjnmqaeeUkREhK666ipVVVVpzJgxevnllxUeHm6Oee211zRt2jRzNdJJkyZp3rx5HXewsDQCYdtROwAAYBfchzDEcY8WNOdblOjw4cMnXMkWp1dVVaW4uDhJhEMAAAKB17gdh2sIAYsiyLQdtQMAAHZBIAQsilDjHw0NDYGeAgAAQLshEAJAM43DNMEaAABYGYEQsCiCjH/QIQQAAFZGIAQsikDYdo1rRyAEAABWRiAELIpA6B8EQgAAYGUEQgBohg4hAACwCwIhYFF0CP2DQAgAAKyMQAhYFIGw7egQAgAAuyAQAsApEAgBAICVEQgBi6JD2HZ0CAEAgF0QCAGLIhD6B4EQAABYGYEQsCgCYdvRIQQAAHZBIASAUyAQAgAAKyMQAhZFh9A/CIQAAMDKCISARREI245TRgEAgF0QCAHgFAiEAADAygiEgEXRIWw7OoQAAMAuCISAhTQOMgRC/yAQAgAAKyMQAkAzdAgBAIBdEAgBC6FD6H8EQgAAYGUEQsCiCIRtR4cQAADYBYEQsBA6hP5HHQEAgJURCAGgGTqEAADALgiEgIXQIfS/+vr6QE8BAACg3RAIAQshEPoHdQQAAHZBIASAUyAQAgAAKyMQAhZCZ8s/qCMAALALAiFgUQQZ/6COAADAygiEgIXQ2fIP6ggAAOyCQAgAzRAIAQCAXRAIAQshyPgfdQQAAFZGIAQshEDoH9QRAADYBYEQAJppHAIbGhoCOBMAAID2RSAELITOln9QRwAAYBcEQsBCCDL+Rx0BAICVEQgBCyG8+AfBGgAA2AWBELAogkzbEQgBAIBdEAgBCyHI+B91BAAAVkYgBCyEQOgf1BEAANgFgRCwEMKLf3DbCQAAYBcEQsBC6Gz5B3UEAAB2QSAELIQg43/UEQAAWBmB8DTmzJkjh8Oh3Nxcc5thGJo5c6bcbrdiY2M1evRobdu2rcnXeb1eTZ06VWlpaYqPj9ekSZN04MCBJmOKi4uVk5Mjp9Mpp9OpnJwclZSUdMBRATgVgjUAALALAuEpbNy4US+88IJ++MMfNtn++OOPa+7cuZo3b542btwol8ulsWPHqqyszByTm5urJUuWaPHixVq9erXKy8s1ceJE1dfXm2MmT56szZs3a+nSpVq6dKk2b96snJycDjs+WA9Bxj+oIwAAsAsC4UmUl5fr2muv1YIFC5ScnGxuNwxDTz/9tB544AH9/Oc/V3Z2tl555RVVVlZq0aJFkiSPx6MXX3xRTz75pC6++GINHjxYCxcu1JYtW/Thhx9KkrZv366lS5fq//7f/6vhw4dr+PDhWrBggd555x3t3LkzIMeM0EeQ8T/qCAAArIxAeBK33XabLrvsMl188cVNtufn56uwsFDjxo0zt0VHR2vUqFFas2aNJCkvL0+1tbVNxrjdbmVnZ5tj1q5dK6fTqWHDhpljzj//fDmdTnMM8H0RCP2DOgIAALuICPQEgtHixYu1adMmbdy4scVjhYWFkqSMjIwm2zMyMrRv3z5zTFRUVJPOom+M7+sLCwuVnp7e4vnT09PNMSfi9Xrl9XrNz0tLS1t5VLADwot/EAgBAIBd0CFspqCgQHfccYcWLlyomJiYk45zOBxNPjcMo8W25pqPOdH40z3PnDlzzEVonE6nMjMzT7lP2BdBpu24DyEAALALAmEzeXl5Kioq0pAhQxQREaGIiAitXLlSzzzzjCIiIszOYPMuXlFRkfmYy+VSTU2NiouLTznm8OHDLfZ/5MiRFt3Hxu6//355PB7zo6Cg4IyOF9ZCZ8v/qCMAALAyAmEzY8aM0ZYtW7R582bzY+jQobr22mu1efNm9erVSy6XS8uXLze/pqamRitXrtSIESMkSUOGDFFkZGSTMYcOHdLWrVvNMcOHD5fH49GGDRvMMevXr5fH4zHHnEh0dLQSExObfAA+hBf/IFgDAAC74BrCZhISEpSdnd1kW3x8vFJTU83tubm5mj17tvr06aM+ffpo9uzZiouL0+TJkyVJTqdTN9xwg6ZPn67U1FSlpKTo7rvv1qBBg8xFavr3768JEyZoypQpev755yVJN954oyZOnKi+fft24BHDSggy/kEdAQCAXRAI22DGjBmqqqrSrbfequLiYg0bNkzLli1TQkKCOeapp55SRESErrrqKlVVVWnMmDF6+eWXFR4ebo557bXXNG3aNHM10kmTJmnevHkdfjywDoKM/1FHAABgZQ6DVzshrbS0VE6nUx6Ph9NHoYMHD6pr166SpE8++UQjR44M8IxC0/bt2zVgwABJ0ptvvqmf/vSnAZ4RAAD2wmvcjsM1hIBF8V5P29FpBQAAdkEgBCyEIOMf3HYCAADYBYEQsBACof9RRwAAYGUEQsBCCC/+QbAGAAB2QSAELIQg4x/UEQAA2AWBELAQgoz/UUcAAGBlBEIAaIZgDQAA7IJACFgIQcY/qCMAALALAiFgIQQZ/6COAADALgiEgIUQZPyP+xACAAArIxACFkII9A+CNQAAsAsCIWAhBBn/oI4AAMAuCISARRFk/IM6AgAAKyMQAhZCePEPOoQAAMAuCISAhRBk/IM6AgAAuyAQAhZCkPEP6ggAAOyCQAhYCEHG/7jtBAAAsDICIWAhhED/IFgDAAC7IBACFkWQaTsCIQAAsAsCIWAhBBn/C7Y6VlVVqaioKNDTAAAAFkEgBCwk2MJLqArmYH355ZcrIyMj0NMAAAAWQSAELCSYg0woCeY6rlixItBTAAAAFkIgBCwkmINMKKGOAADALgiEAHAKBEIAAGBlBELAQuhs+Ufj2gXbfQgdDkegpwAAACyEQAhYCIHQP4K5jmFh/NoGAAD+wysLwEKCOciEKuoIAACsjEAIWAjhxT8I1gAAwC4IhICFEGT8I5jryDWEAADAnwiEgEUFW5AJJcEcCAEAAPyJQAhYCOHF/4K1psE6LwAAEFoIhICF0Nnyj1CoY7DdDgMAAIQmAiFgIaEQZEJBMN+H0Ke+vj7QUwAAABZAIAQshEDof8Fax2ANqgAAILQQCAELCdbwEmqCOVj7VhklEAIAAH8gEAIWFWxBJpQEcyD0IRACAAB/IBACFhIKQSYUhEIduYYQAAD4A4EQsJBgDS+hLFhrSocQAAD4A4EQsJBQ6GyFglCoI4EQAAD4A4EQsJBQCDKhIJhvO8GiMgAAwJ8IhICFEAL9L1hrSiAEAAD+QCAELCpYg0woCIVOK4vKAAAAfyAQAhYSCkEmFIRCHekQAgAAfyAQAhYSCkEmFPhq53A4graOBEIAAOAPBMJm5syZo/POO08JCQlKT0/XFVdcoZ07dzYZYxiGZs6cKbfbrdjYWI0ePVrbtm1rMsbr9Wrq1KlKS0tTfHy8Jk2apAMHDjQZU1xcrJycHDmdTjmdTuXk5KikpKS9DxEWFqzhJVQFYyBkURkAAOBPBMJmVq5cqdtuu03r1q3T8uXLVVdXp3HjxqmiosIc8/jjj2vu3LmaN2+eNm7cKJfLpbFjx6qsrMwck5ubqyVLlmjx4sVavXq1ysvLNXHixCbX/UyePFmbN2/W0qVLtXTpUm3evFk5OTkderywFjqE/uGrXVhYWNDVkUAIAAD8KSLQEwg2S5cubfL5Sy+9pPT0dOXl5enCCy+UYRh6+umn9cADD+jnP/+5JOmVV15RRkaGFi1apJtuukkej0cvvviiXn31VV188cWSpIULFyozM1Mffvihxo8fr+3bt2vp0qVat26dhg0bJklasGCBhg8frp07d6pv374de+CwBAKhf4RCIGRRGQAA4A90CE/D4/FIklJSUiRJ+fn5Kiws1Lhx48wx0dHRGjVqlNasWSNJysvLU21tbZMxbrdb2dnZ5pi1a9fK6XSaYVCSzj//fDmdTnPMiXi9XpWWljb5ANA+HA5H0HbignVeAAAgtBAIT8EwDN11110aOXKksrOzJUmFhYWSpIyMjCZjMzIyzMcKCwsVFRWl5OTkU45JT09vsc/09HRzzInMmTPHvObQ6XQqMzOz7QcIy6FD6B/B3CH0IRACAAB/IBCewu23364vvvhCf//731s85jtty8cwjBbbmms+5kTjT/c8999/vzwej/lRUFBwusOAjRAI/SOYAyHXEAIAAH8iEJ7E1KlT9dZbb+njjz9Wt27dzO0ul0uSWnTxioqKzK6hy+VSTU2NiouLTznm8OHDLfZ75MiRFt3HxqKjo5WYmNjkA/AhEPqHr3bh4eFBV0cCIQAA8CcCYTOGYej222/XG2+8oY8++khZWVlNHs/KypLL5dLy5cvNbTU1NVq5cqVGjBghSRoyZIgiIyObjDl06JC2bt1qjhk+fLg8Ho82bNhgjlm/fr08Ho85Bvi+gi28hLpgvu0Ei8oAAAB/YJXRZm677TYtWrRI//rXv5SQkGB2Ap1Op2JjY+VwOJSbm6vZs2erT58+6tOnj2bPnq24uDhNnjzZHHvDDTdo+vTpSk1NVUpKiu6++24NGjTIXHW0f//+mjBhgqZMmaLnn39eknTjjTdq4sSJrDAKvwi2IBNKgvmUUR86hAAAwB8IhM3Mnz9fkjR69Ogm21966SX9+te/liTNmDFDVVVVuvXWW1VcXKxhw4Zp2bJlSkhIMMc/9dRTioiI0FVXXaWqqiqNGTNGL7/8ssLDw80xr732mqZNm2auRjpp0iTNmzevfQ8QlsYpo/4RzIGQU0YBAIA/EQibac2LP4fDoZkzZ2rmzJknHRMTE6Nnn31Wzz777EnHpKSkaOHChW2ZJnBCwRZeQh23nQAAAFbHNYSAhdAh9A86hAAAwC4IhICFEAj9IxQCIYvKAAAAfyAQAhYSbOElVAXzbSd86BACAAB/IBACFhWsQSaUBPNtJwiEAADAHwiEgIVwyqh/BPMpoz4EQgAA4A8EQsBCCIT+EcyB0Nch3L59e4BngmBXVFQU6CkAAEIAgRCwkGALL6EuGE8Z9bntttsCPQUEsXfffVcZGRnaunVroKcCAAhyBELAQhqfRhisQSYUNO4QBtupmb4OIXAqmzdvliTl5+cHdiIAgKBHIAQspPGtCAiEbRfMp4wCAAD4E4EQsBDuTecfwXzbidra2kBPASHA10kOtu9fAEDwIRACFsIpo/4VjNcQVlVVBXoKCAEEQgBAaxEIAQvhlFH/COZTRnv27ClJuvzyywM7EQAAYAkEQsBCfIEwGDtbocRXu2CsY1jY8V/bbrc7wDNBMKNDCABoLQIhYCG+U0bDw8MDPJPQ5qtjRERE0L2g9s2nrq4uwDNBMCMQAgBai0AIWIivQxiMi6GEEl8dgzEQ+sIqi8ugNYLt+xcAEHwIhICF+MJCMF77FkoaB8Jguw+hHTqEK1assPTxdQTuVwkAaC0CIWAh9fX1CgsL48XgGQqFDqFVA9POnTt18cUX68knnwz0VEKa73cAt6IBAJwOgRCwkPr6evP6wWALMqGkcbAOtjpavUNYXl4uSSooKAjwTKyBU4sBAKdDIAQspKGhQeHh4UEZZEKJL1gHYx2tfg2hbxXVYDtVN9T4OoRWfeMAAOA/BELAQjhl1D9CIRBa9YU+gdC/rPrGAQDAfyICPQEA/sMpo/4RzIHQ6qeM+t7MIBCeGat3kgEA/kOHELAQ3ymjEoHwTARzILT6C33f8REIz4zvDQOrfp8AAPyHQAhYSDAvhhJKfIEwLCws6IKJ7//VqqtH+oJMoOpuGIYlfnYIhACA1iIQAhbSeFEZtF0odAitespooAPhT37yEyUkJARk3/7kq2NNTU2AZwIACHZcQwhYiK9DKHHK6JkI5tVarX4NYaAD4cqVKwOyX3+jQwgAaC06hICFBHNnK5QE86m3dAjRGnQIAQCtRSAELIRTRv0jmIO1LyhZ9RpCX0eLQHhm6BACAFqLQAhYCKeM+kcwB0JOGUVr+Or42GOPacWKFQGeDQAgmBEIAQsJ5iATSmpraxUZGRmUdWxoaJDD4SAQ4pQaf3/84x//COBMAADBjkAIWEgwL4YSSqqrqxUTE6OwsLCgq6NhGIqMjLTsKaMEQv9ofKqo1+sN4EyAk5s7d65cLlegpwHYHquMAhbSeDEUtJ3X61VMTIwcDkfQBZP6+npFRUXRIcQpNQ6ELCyDYDV9+vRATwGA6BACluLrEErBcQ3hokWLtGrVqkBP43vzdQiDsdNaW1ur2NhYAiFOqXEg7MhafvXVV3riiSc6bH9W9uqrr6q4uDjQ02hXvjcv6WIDgUUgBCyktrZWERERQRNkrr32Wo0aNSrQ0/jegjUQGoahuro6AmEHCPT+z1SgAuHkyZM1Y8YMy35/dhSPx6PrrrtON9xwQ6Cn0q58i6CVlJQEdiKAzREIAQtpHGTQdr46hoeHB9W1er4X2bGxsUE1L3/yHePbb7+to0ePBmweod6xCFQg9L2BcuzYsQ7bpxX56ldYWBjgmbQv3xktBEIgsAiEgIX4gowUHKeMhipfHWNjY1VdXR3o6Zh814LZoUMoSR9//HHA5lFWVhawfftD4+sGOzIQ+n7/BDLMW0F5ebkk695eprlQ/3kDQh2BELAQr9er6OjooDvVMdT4AmFcXJwqKioCPR2Tr2sVHx9v2ReKjY8rMjIyYPPo379/wPbtD4HqEPoCoRU6hEePHtW3334bkH37fu80/n+0It/fqWD6PQvYEauMAhbCKaP+0TgQVlZWBno6Jl8gjIuLs0Ug7OgXw41Pww31QNO4dlVVVR22X981Yb4OVyhLS0tTQkKCSktLO3zfdgmEvjcrrPD9AoQyOoSAhfg6hAkJCW2+JsNf3YRQ7lAGeyCMj4+3/DWEUscGmUDsrz3V1NRo2LBhcjqdHdp98b0ZZZWOT6BOZbRjILT6sQLBjEAIWEhFRYXi4+OVlpbWpuXKlyxZovj4eL+8Ix5M1959X75gHcyB0KodwkB1tgKxv/ZUWVmpIUOG6Je//GVAwtnOnTst+6ZFR/D9n1n151w6/qah743DlStXKioqSv/4xz8CPCvAngiEgIWUlpYqISFBiYmJbQp1ixcvVnV1tfbt23fGc2n84jrUboxdXl6uuLg4xcfHB1Wno3EgNAzDki+4q6qqFB8fb/67o/dtFZWVlQH5HvYF+t///veaNWtWh+23PQWic+X7P/N4PIqPj1d+fn6Hz6G9NX6z7YsvvpAkvfDCC4GaDmBrBELAQsrKypSYmCin0ymPx/O9vz4qKkpS65c6f+GFF7R69eoTPtb4xXUgrsFpq/r6ehUWFsrtdishIUHl5eVBE7x8gTApKUmSgqp76S9VVVVyuVxKSUkJeCAM5dOeAxUIG39Prl+/vsP262+NO3OHDx/u8P37/s+OHj2qyspKjRkzpsPn0N4an8XiOzU3IyMjUNNp4ciRIxozZowKCgoCPRWg3REIAQspKytTQkKCUlNTdeTIke/99b5A2NoXQDfddJN+/OMfn/Cxxi8M2xJOA+XQoUOqr69XZmamUlNTZRiGPv/880BPS9L/1rFr166Sjv9/b9iwQd98800gp+VXviATGxvb4YHQ9z171113Nfk8FFVUVAQkEDbeV1xcXIft198a/87KzMzUxo0bO3T/za9dHDduXIfuvyMcOHDA/PeePXskBdd9F9euXauPPvpIjz32WKCnArQ7AiHazDCMkF+Jz0oMwzBPGe3SpYsOHTr0vZ/Dd91fa/4on657EqodQt+LlG7duik1NVWSNGTIkKDoEhYVFUmSevfuLen4i8Zhw4apZ8+eHbL/VatW6ZZbbmnXfQRDIExPT5cUuvdGa2hokMfjkdPpDGggjI2N7bD9+lvzRbl+9KMfdej+m9/H0Yr3ddy/f7+k47X1/ezt3r07kFNqwnfLESu94QacDIEQbfbwww8rNTVVDodDDodDKSkp+t3vfqc777yzzStcou327Nmj2tpa9erVS263W0eOHDFPMWwtX8BvTSA83TLhjV8YhlKH8F//+pek4104p9Npbn///fcDNSWT7xS8Ll26SPrfwNJRC0/87Gc/03PPPef3a0Lz8/PlcDj04YcfyuPxKCEhISCB8ODBg5Kks846S1LoBsLS0lI1NDQoJSVF8fHxqq2t7bDr4Br/3IfSG0HNNf8b9tOf/rRD93/06FF16tTJ/Nwf13UHm4KCAnXq1ElZWVlNtgXLgmQffvihJGnLli0yDEO33HKL1q1bF+BZAe2DQBgE/vrXvyorK0sxMTEaMmSIPvnkk0BPqVX69evX5PPi4mK9+OKLevrpp5WcnGwGRd/HjTfeqJdeeknz5s3TX/7yFy1atEj79+9XYWFhSF+rEyx+//vfS5KGDx+uPn36yDAM8zSc1vJd09GaQNj4lNQTdc8an3YaSp3kRx99VJKUmpqqc845x9x++eWXB2hG/+uZZ56RJCUmJkpqeo2W79329hQdHS1J+vjjj/36vO+9954k6ZZbbtG3336rtLQ0JSYmat68ed/7TY0zceDAAcXFxal79+6STh0IDcMI2sWSfF2Wrl27aufOnZKkvLy8dt9vfX29PB6P+XPTEd+T7aV5IOzo++QtXLhQgwcPNj+3YiBcvHixIiMjNXDgQElSQkKCDMPQu+++G+CZHef7mfnqq69011136bnnntPw4cMDPCugfRAIA+z1119Xbm6uHnjgAX322Wf68Y9/rEsuuSQk/pBOnjxZBw4c0G9/+9tWjV+wYIF++9vfaurUqbr99tt17bXXqkePHurSpYvCwsJaBMgTfaSnp+uuu+7SggULtHDhQu3atUv//Oc/VVhYqDvvvFNvv/22li1bpt27d+vIkSMqLy8P+bBZW1t72k5JYWGh/vnPf0o6vuBI//79JUlTp079Xp0BX4ekcSDctWvXCQOf7/RF6fi7useOHWtymmpBQYEiIyMVGRlpPm9H2rlzp/bu3duqseXl5dq1a5fOPvtsSdLPf/5z835qvu+fsLCwgH4vNX7X3PeO+r///W9z2xtvvCHpeKhv3JlpaGjw27x9NZkwYYJfns/HF6z27Nmjf//738rMzNSnn34qSRo1apRf93UqW7dulcvlMgP3qU7TmzlzpqKjo4PyzQ7f6Y3Z2dn69a9/LUlavnz5CcdWVVX5LeDv2LFDVVVVmjt3rmbNmqWvv/5ahmFo586dWrJkiV/20VE2bdpkvgEi+f9NkFPx1a1r165atWqVHnnkERUVFZ30jTqPxyOHw6Hbb79dtbW1euSRR/Rf//VfHTbftvjmm2/0n//8R8XFxXrwwQf1wQcfmIu3PP/88wGe3fEQuGvXLvP0/Kefftp8zNc5PJXRo0fL4XCouLhYDQ0NQbVaNXAiDiPUXy2HuGHDhuncc8/V/PnzzW39+/fXFVdcoTlz5pz260tLS80VJX0vYgKtoaFB+/btk9PpVH5+viZPnqyhQ4dq0aJFgZ5aUOnVq5e+/vprxcbGKiEhoUnA6tKliwYPHqxPPvnE7FKkpqae9jqSp59+WnfccYckyeVytVgc5vrrr9eWLVvk9Xp10UUX6V//+peOHTt2wne/Bw4cqG3btpmfJycnq7i4WD/5yU+UlZWl//7v/zYf69Onj9mVGDlypLnyaGZmpvlHPj09Xd27d5fb7VZDQ4Peeecd3Xfffbrgggv0zjvvqHPnzurSpYu2bdumn/70p9q1a5f279+v0tJSRUVFqba2VhdeeKEKCgr06aef6q233tL06dPVqVMn/fGPfzTnkZCQoL59++rvf/+7pOPX/w0dOlR///vfW30KW0lJSZPTRefNm6epU6c2GRMVFaWamhpFR0ebXayoqCglJyerV69e+tGPfqS6ujp1795dvXv31o4dOxQREaHCwkJ169ZN06dPlyRddNFFamho0FlnnaWNGzfqs88+kyQNHjzY/Hdj//M//6Nf/OIXOv/8880OYVhYmHmDZ5/x48crPj7eDIpXXnmlhg0bpjFjxqi6ulr19fVKSkpSQkKCvv32W8XExKiyslKLFi3S9u3b9eCDDyoqKkrx8fGqrq7Whx9+qPvvv998/iVLlqhHjx7q27ev9u3bpx/84AcyDEP79+9XWVmZPB6PunbtquTkZFVUVGjbtm2aOHGipOOLEfXp00eHDx/WRRddpEsuuaTJ3NevX689e/bo2muvlXT8es65c+eat9rIzs5WfHy8ysrKVFtbq4qKCg0cOFCbNm2S2+3W3r17NWvWLKWnp+vOO+9Uv379VFhYqD179igtLU2HDh3Spk2bVFhYqJEjR+pPf/pTkwUuamtrFRkZKanlz0F+fr42bdqkK6+80ty2Zs0a5efn69JLL1VpaakKCwvlcrnMn5mYmBgVFRWpuLjYfDE5YcIEhYWFyev16ssvv1RERITGjh0rj8ej5ORkGYYhr9er2NhYRUZGqrS0VK+//rr69u2rgoICuVwuVVRU6NFHH9XQoUMVExOjv/3tb03qaBiGGhoaFB4e3mT7n//8Z1VWViopKUnTpk0zt19zzTUaNWqUevXqpd69e+uTTz5R586dlZSUpE2bNqlz586KiopSaWmpXnrpJTkcDkVFRWnQoEFauHCh+fumtLRUq1atMv+/G7vooov03HPPadWqVRowYIDq6+uVmJiosLAw9ezZUwcOHFBycrK2bt0qh8OhMWPGaNasWQoLC9Oll14q6fhiNS6XS1VVVSouLpZhGCopKdGOHTs0dOhQ7dmzR127dlVpaan5/de5c2f169dPPXr00MCBAzV+/HjFxcVpypQpkiS326358+dry5YtevDBByUd/x29detWvf766/rNb34jSZo/f75cLpeeeeYZXXHFFRo6dKh27dqliy66SF9//bX+8Y9/aNSoUerXr5927dqlXr16qby8XNu3b9eqVau0aNEiderUSSNGjNDYsWOVnJysvXv3asuWLebp6j5fffWVevXqpd27d5unMV955ZWKiIhQ165dVVZWptdff/20v9euvvpq/fSnP9X555+vtWvXKiUlRQkJCebvrvj4eO3YsUORkZHKyMjQ1q1bdeGFF6q0tFSbNm1SeHi4unbtqrVr12r27NmSjq8GWlNToylTpqiqqkrXXHONPv/8c0VFRamoqEgul0vl5eX66KOPlJGRoW7duumRRx5p8Qblww8/rJkzZ5qf+9508omMjNSFF16o7OxsHTx4UOecc452796tgoICrV+/Xn/5y19UWFioSy+9VKtXrzYXBHvyySd19dVXq7i4WEeOHFFVVZWSk5M1YMAAJSQkqKysTF6vV//zP/+jgoICs5Pe3JIlS/Szn/1Mklpcj5uUlKSSkhL97Gc/U8+ePfXVV19pw4YNOuuss7Rq1aoTPt8FF1ygH/7wh/r222914YUXat++fYqPj1dSUpLy8vI0YcIElZeX69ChQ/ryyy8VHR2tvXv3avXq1Ro5cqTuuOMOHT58WJs2bVLPnj116NAhlZWVqXfv3urVq5feeOMNXXfddTp27Jg+//xznXPOOerUqZO6dOmitWvXauDAgUpPT9ctt9yiiy66SNnZ2fr8888VExOjXbt2afTo0aqurtbu3bvldDr13nvvafbs2QE7QyYYX+NaFYEwgGpqahQXF6d//OMf5i8cSbrjjju0efNmrVy5ssXXeL3eJqdQlZaWKjMzMzA/LK+/Lr32Wpu/vK6uThEREcf/XV8vT0mJSsvKFBEerpqaGtXX16u8vFyVlZUqLS1VXX294uLiVF9f36GnkYWK9M6ddf7555uf1zc0tOnUm+zsbG3durXV431/XE+kX79+2rFjx/eeQyBNmDBBUd+FAZ+Ghga9EwSnMfX+wQ80YMAAScd/F3ywbJliY2LUxe3W119/3SFzOPvss9tt1dXeP/iBKisrNXToUEnSseLik97WpL2ce+656ta1qz5ZvbrJsvihxCFp7Lhxivmuw3XgwAFtOsGbC+0hPi5OY8aMUW1dnZa+/75C+QVG7969NaB/fxmGobffeadD952QkKCfjB5tfv7ee++pLggWtvKnrJ49NWjQoCbb9u7dqy+2bAnQjJrKHjhQvXr1Ul19vY4dParUtDR9tmmTDrZhwbZQ1iDpx99+ay6y1pEIhB2HQBhABw8eVNeuXfXpp59qxIgR5vbZs2frlVdeOeE7VjNnzjS7IY2FYiAMJoZhqLKyUhGRkcfDgMOhutpaOcLCFB4Wpvr6elVWVSkuNlZh4eEyGhokh0M1Xq9q6+pUV1uryKgoHTt2TCnJyaqqqlJcXJyqvV6VlZaqvKJCCZ06KSEhQQ2GobraWkVERKiyqkrJSUn66quvlJGRocOHD6uuvl49uneX1+uVw+FQUlKSYmJitHfvXjUYhlJTUlRTU6OIiAjFd+qkqKgoRTTrADRW/90pg4Xf/RErLi5WfUODarxelZeXq7yiQrGxserWtav69OmjiIgIlZaVafeuXerSpYtSUlJ06NAhHT58WJWVlTqrb19VV1Wpa9euiomJkWEY2r17t+I7dVJycrI8JSWSw6HY2FglJyUdX/GwtFRHjx7V3vx8GTp+mtqIESPM7/FO8fHat3+/Mrt1kyMsTOXl5QoLC9OxY8fU0NAgZ2Ki6hsaFBkRIUMyH8vq2VPVXq+cTqdiY2JUXFys2ro6lZeX66yzzlKn+Pjjcy8qUnJysiLCwxUVFaUePXuqqrJS0dHRioiMlOOk1WtZy+rqakVFRSnM4VBJSYlq6+pUXFys3bt3KysrS/Hx8cdXpzMMRUZF6UBBgRKdTkWEhysiMlJGQ4NiY2PldDqVnpGhj1asUHhEhBoaGtRQXy93167yer1KcjqVmJio8IgIderUSZ3i41u8e97YwYMHVeLxqEf37goLC9N/8vJUU1Oj7OxsOZ1O/fvjj1VTW6vu3bur8NAh1dXXq6GhQXFxceYKf76uZ1JSkmq+e9Ol8rvTlWOio3X2OecoIz1du/fs0f59+1RRWalOnTopJjpamZmZqqqu1o4dO1p0s6MiI9XQ0KABAwcqyenU9h07dPTbb9VgGOrqdsvpdKpnz57mG0SNNTQ0aNuXX+rbI0dUVl6uyMhIJSUlKSUlRcXHjikyMlIlJSWKiIxU1Xfz8dbUqKKiQllZWTIaGrR//37zZ+fod6d4+t7ISE5KUmlZmfr17asePXo0mYOv69QzK0tpqanav3+/qqqqVHTkiFwZGerXv78OFBTos82bJUk9e/TQ3n37lJyUpPr6enXu3FlFR47IW12tmu86Iw4d764cO3ZMNbW1io6Kkve7U2ZTU1PNn4WY6GhVf/d/4Ha7JcMwX4h2795d+/fvV1JSkpKTk1VfX6/0zp1lGIZSU1MVExPT8nu3vl4VFRUqLStTfV3d8b8ZTqeKi4s1oH9/1dXXa9++fSopLlZVdbXZCfHNLzwsTBkZGYqKjlaS06nCw4cVFxuraq9Xx44elbtrV/Xp3bvJaZaSVFxSosqKCnXt2lXFxcXauXOnkpKTVVRUpJKSEsXExCgiIkLl5eWKj4s7/j0VH6/y7/Y/+JxzVFFZqYqKCh07elSGpPCwMNV8t1BObGysYmNi5PV6VVFZqS4ulw41OrXSlZEhQ1JERIRSUlJUUVGh4uJiFRcXKyoqSj867zxVVlYqNi5Oe/bsUX19vdxutxITE5WSnNzk+/Do0aPavWeP6urqVFJSYv7OdrvdciYm6lBhoXntYUREhOrq6pSamqramhpVVlWpT+/e6ty5s3bv2aPOaWnyeDyqqqpSZVWVMjMzZRiGkpxOpaSknPBnob6+XqWlpWYHPjIyUnHx8YqOjlZaWprq6+sVHh6u8LAwVVRWqrq6WlVVVfr222+1f/9+9T3rLB08dEhpqakqLCyU1+tVVlaWDhcVqaK8XKlpaUro1EmFhw+rqqpKsbGxiouNVWVVlXnZQnJysgYPHqyvv/5aXd1ufbFliyorKtS9Rw/V1NSorq5ONTU1Ki4uVlpamiIjIlRTW6uEhATFx8crJjpaqWlp5u/BE6mrrz/+N7eyUvv27lVcfLy81dU6/F3nMTo6WkWHD+vosWMKCwtTV7dbcfHxKvV4VOLxKDoqSknJyebnLpdL1VVVCv/uzeaY2Fh5q6tlGMbxBfFSU9WjR4/jp4AfParO6ekKDzv5FVVer1elZWVy6PjlFWXl5erSpYtKv1s8LTU1VW63Ww6HQ9XV1TIkxcbEqLCwUF988YUcYWGKi41VTGysjhQVqaa2VrExMaqqrlZWVpY8JSWK+u509E6dOikyIkJd3G59c+CA+vTpo2PfvVEVER6uggMHlNCpk2q/e01hSEro1En7CwrM/7NeWVmKi4uTp7RU4WFhxzv5DodcLpdiY2JU4vEc/5v73eKA5eXlqquvlysj4/jv0d69lbhixUnr0Z4IhB2HQBhAvkC4Zs2aJhcqz5o1S6+++uoJOytB1SEEAAAA2gGBsOO0fAsKHSYtLU3h4eEtLhQvKipSRkbGCb8mOjq6xTuwAAAAANAWrDIaQFFRURoyZEiL1d+WL1/e5BRSAAAAAGgPdAgD7K677lJOTo6GDh2q4cOH64UXXtD+/ft18803B3pqAAAAACyOQBhgV199tY4ePapHHnlEhw4dUnZ2tt577z316NEj0FMDAAAAYHEsKhPiuOAWAAAAVsNr3I7DNYQAAAAAYFMEQgAAAACwKQIhAAAAANgUgRAAAAAAbIpACAAAAAA2RSAEAAAAAJviPoQhznfXkNLS0gDPBAAAAPAP32tb7pDX/giEIa6srEySlJmZGeCZAAAAAP5VVlYmp9MZ6GlYGjemD3ENDQ06ePCgEhIS5HA4TjimtLRUmZmZKigo4Mae7YD6th9q276ob/uhtu2L+rYfatu+qG/rGYahsrIyud1uhYVxlVt7okMY4sLCwtStW7dWjU1MTOSXTzuivu2H2rYv6tt+qG37or7th9q2L+rbOnQGOwZxGwAAAABsikAIAAAAADZFILSB6OhoPfzww4qOjg70VCyJ+rYfatu+qG/7obbti/q2H2rbvqgvghGLygAAAACATdEhBAAAAACbIhACAAAAgE0RCAEAAADApgiEAAAAAGBTBEIb+Otf/6qsrCzFxMRoyJAh+uSTTwI9pZAzZ84cnXfeeUpISFB6erquuOIK7dy5s8kYwzA0c+ZMud1uxcbGavTo0dq2bVuAZhy65syZI4fDodzcXHMbtT0z33zzjX71q18pNTVVcXFxOuecc5SXl2c+Tn3bpq6uTg8++KCysrIUGxurXr166ZFHHlFDQ4M5htq23qpVq3T55ZfL7XbL4XDozTffbPJ4a2rp9Xo1depUpaWlKT4+XpMmTdKBAwc68CiC06lqW1tbq3vvvVeDBg1SfHy83G63rrvuOh08eLDJc1Dbkzvd925jN910kxwOh55++ukm26kvAolAaHGvv/66cnNz9cADD+izzz7Tj3/8Y11yySXav39/oKcWUlauXKnbbrtN69at0/Lly1VXV6dx48apoqLCHPP4449r7ty5mjdvnjZu3CiXy6WxY8eqrKwsgDMPLRs3btQLL7ygH/7wh022U9u2Ky4u1gUXXKDIyEi9//77+vLLL/Xkk08qKSnJHEN92+axxx7Tc889p3nz5mn79u16/PHH9cQTT+jZZ581x1Db1quoqNDZZ5+tefPmnfDx1tQyNzdXS5Ys0eLFi7V69WqVl5dr4sSJqq+v76jDCEqnqm1lZaU2bdqkhx56SJs2bdIbb7yhXbt2adKkSU3GUduTO933rs+bb76p9evXy+12t3iM+iKgDFjaj370I+Pmm29usq1fv37GfffdF6AZWUNRUZEhyVi5cqVhGIbR0NBguFwu49FHHzXHVFdXG06n03juuecCNc2QUlZWZvTp08dYvny5MWrUKOOOO+4wDIPanql7773XGDly5Ekfp75td9lllxm//e1vm2z7+c9/bvzqV78yDIPanglJxpIlS8zPW1PLkpISIzIy0li8eLE55ptvvjHCwsKMpUuXdtjcg13z2p7Ihg0bDEnGvn37DMOgtt/Hyep74MABo2vXrsbWrVuNHj16GE899ZT5GPVFoNEhtLCamhrl5eVp3LhxTbaPGzdOa9asCdCsrMHj8UiSUlJSJEn5+fkqLCxsUuvo6GiNGjWKWrfSbbfdpssuu0wXX3xxk+3U9sy89dZbGjp0qH7xi18oPT1dgwcP1oIFC8zHqW/bjRw5UitWrNCuXbskSZ9//rlWr16tSy+9VBK19afW1DIvL0+1tbVNxrjdbmVnZ1Pv78nj8cjhcJhnElDbM9PQ0KCcnBzdc889GjhwYIvHqS8CLSLQE0D7+fbbb1VfX6+MjIwm2zMyMlRYWBigWYU+wzB01113aeTIkcrOzpYks54nqvW+ffs6fI6hZvHixdq0aZM2btzY4jFqe2a+/vprzZ8/X3fddZd+//vfa8OGDZo2bZqio6N13XXXUd8zcO+998rj8ahfv34KDw9XfX29Zs2apV/+8peS+N71p9bUsrCwUFFRUUpOTm4xhr95rVddXa377rtPkydPVmJioiRqe6Yee+wxRUREaNq0aSd8nPoi0AiENuBwOJp8bhhGi21ovdtvv11ffPGFVq9e3eIxav39FRQU6I477tCyZcsUExNz0nHUtm0aGho0dOhQzZ49W5I0ePBgbdu2TfPnz9d1111njqO+39/rr7+uhQsXatGiRRo4cKA2b96s3Nxcud1uXX/99eY4aus/bakl9W692tpaXXPNNWpoaNBf//rX046ntqeXl5en//qv/9KmTZu+d62oLzoKp4xaWFpamsLDw1u8u1RUVNTiXVa0ztSpU/XWW2/p448/Vrdu3cztLpdLkqh1G+Tl5amoqEhDhgxRRESEIiIitHLlSj3zzDOKiIgw60dt26ZLly4aMGBAk239+/c3F5bie7ft7rnnHt1333265pprNGjQIOXk5OjOO+/UnDlzJFFbf2pNLV0ul2pqalRcXHzSMTi52tpaXXXVVcrPz9fy5cvN7qBEbc/EJ598oqKiInXv3t38G7dv3z5Nnz5dPXv2lER9EXgEQguLiorSkCFDtHz58ibbly9frhEjRgRoVqHJMAzdfvvteuONN/TRRx8pKyuryeNZWVlyuVxNal1TU6OVK1dS69MYM2aMtmzZos2bN5sfQ4cO1bXXXqvNmzerV69e1PYMXHDBBS1ukbJr1y716NFDEt+7Z6KyslJhYU3/jIaHh5u3naC2/tOaWg4ZMkSRkZFNxhw6dEhbt26l3qfhC4O7d+/Whx9+qNTU1CaPU9u2y8nJ0RdffNHkb5zb7dY999yjDz74QBL1RRAI0GI26CCLFy82IiMjjRdffNH48ssvjdzcXCM+Pt7Yu3dvoKcWUm655RbD6XQa//73v41Dhw6ZH5WVleaYRx991HA6ncYbb7xhbNmyxfjlL39pdOnSxSgtLQ3gzENT41VGDYPanokNGzYYERERxqxZs4zdu3cbr732mhEXF2csXLjQHEN92+b66683unbtarzzzjtGfn6+8cYbbxhpaWnGjBkzzDHUtvXKysqMzz77zPjss88MScbcuXONzz77zFzpsjW1vPnmm41u3boZH374obFp0ybjoosuMs4++2yjrq4uUIcVFE5V29raWmPSpElGt27djM2bNzf5G+f1es3noLYnd7rv3eaarzJqGNQXgUUgtIG//OUvRo8ePYyoqCjj3HPPNW+VgNaTdMKPl156yRzT0NBgPPzww4bL5TKio6ONCy+80NiyZUvgJh3CmgdCantm3n77bSM7O9uIjo42+vXrZ7zwwgtNHqe+bVNaWmrccccdRvfu3Y2YmBijV69exgMPPNDkRTS1bb2PP/74hL9nr7/+esMwWlfLqqoq4/bbbzdSUlKM2NhYY+LEicb+/fsDcDTB5VS1zc/PP+nfuI8//th8Dmp7cqf73m3uRIGQ+iKQHIZhGB3RiQQAAAAABBeuIQQAAAAAmyIQAgAAAIBNEQgBAAAAwKYIhAAAAABgUwRCAAAAALApAiEAAAAA2BSBEAAAAABsikAIAAAAADZFIAQAAAAAmyIQAgAAAIBNEQgBAAAAwKYIhAAAAABgUwRCAAAAALApAiEAAAAA2BSBEAAAAABsikAIAAAAADZFIAQAAAAAmyIQAgAAAIBNEQgBAAAAwKYIhAAAAABgUwRCAAAAALApAiEAAAAA2NT/B8Qqy6A2ltayAAAAAElFTkSuQmCC", "text/html": [ "\n", "
\n", "
\n", " Figure\n", "
\n", " \n", "
\n", " " ], "text/plain": [ "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = PowderPattern()\n", "\n", "p.ImportPowderPattern2ThetaObs(data_file, 4) # skip first 4 lines header\n", "# Copper K-alpha1+alpha2. Use \"Cua1\" for Cu-alpha1 only\n", "p.SetWavelength(\"Cu\")\n", "print(p.GetWavelength())\n", "\n", "p.plot(hkl=True)\n" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "### Add crystalline phases\n", "We assume all structures are known.\n", "\n", "This will update the above plot, though the scales will be incorrect." ] }, { "cell_type": "code", "execution_count": 4, "id": "7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "UnitCell : Aluminium oxide - $-alpha(R -3 c:H)\n", " Cell dimensions : 4.76050 4.76050 12.99560 90.00000 90.00000 120.00000\n", "List of scattering components (atoms): 2\n", "Al1 at : 0.0000 0.0000 0.3522, Occup=1.0000 * 0.3333, ScattPow:Al1 , Biso= 0.0000\n", "O1 at : 0.3067 0.0000 0.2500, Occup=1.0000 * 0.5000, ScattPow:O1 , Biso= 0.0000\n", "\n", "Occupancy = occ * dyn, where:\n", " - occ is the 'real' occupancy\n", " - dyn is the dynamical occupancy correction, indicating either\n", " an atom on a special position, or several identical atoms \n", " overlapping (dyn=0.5 -> atom on a symetry plane / 2fold axis..\n", " -> OR 2 atoms strictly overlapping)\n", "\n", " Total number of components (atoms) in one unit cell : 30\n", " Chemical formula: Al3+0.33 O2-0.5\n", " Weight: 16.9935 g/mol \n", "\n", "UnitCell : Zincite(P 63 m c)\n", " Cell dimensions : 3.24950 3.24950 5.20690 90.00000 90.00000 120.00000\n", "List of scattering components (atoms): 2\n", "Zn at : 0.3333 0.6667 0.0000, Occup=1.0000 * 0.1667, ScattPow:Zn , Biso= 0.0000\n", "O at : 0.3333 0.6667 0.3450, Occup=1.0000 * 0.1667, ScattPow:O , Biso= 0.0000\n", "\n", "Occupancy = occ * dyn, where:\n", " - occ is the 'real' occupancy\n", " - dyn is the dynamical occupancy correction, indicating either\n", " an atom on a special position, or several identical atoms \n", " overlapping (dyn=0.5 -> atom on a symetry plane / 2fold axis..\n", " -> OR 2 atoms strictly overlapping)\n", "\n", " Total number of components (atoms) in one unit cell : 4\n", " Chemical formula: O0.17 Zn0.17\n", " Weight: 13.5652 g/mol \n", "\n", "UnitCell : Calcium fluoride(F m -3 m)\n", " Cell dimensions : 5.46250 5.46250 5.46250 90.00000 90.00000 90.00000\n", "List of scattering components (atoms): 2\n", "Ca1 at : 0.0000 0.0000 0.0000, Occup=1.0000 * 0.0208, ScattPow:Ca1 , Biso= 0.0000\n", "F1 at : 0.2500 0.2500 0.2500, Occup=1.0000 * 0.0417, ScattPow:F1 , Biso= 0.0000\n", "\n", "Occupancy = occ * dyn, where:\n", " - occ is the 'real' occupancy\n", " - dyn is the dynamical occupancy correction, indicating either\n", " an atom on a special position, or several identical atoms \n", " overlapping (dyn=0.5 -> atom on a symetry plane / 2fold axis..\n", " -> OR 2 atoms strictly overlapping)\n", "\n", " Total number of components (atoms) in one unit cell : 12\n", " Chemical formula: Ca2+0.021 F1-0.042\n", " Weight: 1.62654 g/mol \n", "\n", "========================== WARNING =========================\n", " In ScatteringPowerAtom::GetTemperatureFactor():\n", " Anisotropic Displacement Parameters are not currently properly handled\n", " for Debye-Waller calculations (no symmetry handling for ADPs).\n", " =>The Debye-Waller calculations will instead use only isotropic DPs\n", "\n" ] } ], "source": [ "for cod_id in cod_phases:\n", " c = CreateCrystalFromCIF(\"http://crystallography.net/cod/%d.cif\" % cod_id)\n", " print(c,\"\\n\")\n", " p.AddPowderPatternDiffraction(c)\n", "\n", "p.FitScaleFactorForRw()\n", "p.UpdateDisplay()\n" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "### Add an automatic background\n", "This uses a Bayesian estimation of the background, which should be\n", "good enough if there is a good separation of the peaks" ] }, { "cell_type": "code", "execution_count": 5, "id": "9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "No background, adding one automatically\n" ] } ], "source": [ "# Add background if necessary\n", "need_background = True\n", "for i in range(p.GetNbPowderPatternComponent()):\n", " if isinstance(p.GetPowderPatternComponent(i), PowderPatternBackground):\n", " need_background = False\n", " break\n", "if need_background:\n", " print(\"No background, adding one automatically\")\n", " x = p.GetPowderPatternX()\n", " bx = np.linspace(x.min(), x.max(), 30) # Number of interpolation points\n", " by = np.zeros(bx.shape)\n", " b = p.AddPowderPatternBackground()\n", " b.SetInterpPoints(bx, by)\n", " # b.Print()\n", " b.UnFixAllPar()\n", " b.OptimizeBayesianBackground()\n", "p.UpdateDisplay()\n" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "### Fit profile, step 1\n", "Conservative fit, starting with a fixed width (W=1e-5) \n", "\n", "Note that as we know the crystalline structures, we don't need to perform a Le Bail fit" ] }, { "cell_type": "code", "execution_count": 6, "id": "11", "metadata": {}, "outputs": [], "source": [ "# Multiple phases, so can't use quick_fit_profile\n", "#\n", "# NOTE: we don't use this function as the phases are known, but that's the way\n", "# to do it with multiple crystalline phases\n", "def do_lebail(init=False):\n", " \"\"\"\n", " This performs a Le Bail fit by looping over all phases,\n", " one at a time. Le Bail is disabled on output\n", " \"\"\"\n", " for i in range(20):\n", " for i in range(p.GetNbPowderPatternComponent()):\n", " pdiff = p.GetPowderPatternComponent(i)\n", " if not isinstance(pdiff, PowderPatternDiffraction):\n", " continue\n", " if i==0 or init:\n", " pdiff.SetExtractionMode(True, True)\n", " else:\n", " pdiff.SetExtractionMode(True, False)\n", " pdiff.ExtractLeBail(1)\n", " pdiff.SetExtractionMode(False, False)\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "12", "metadata": {}, "outputs": [], "source": [ "for i in range(p.GetNbPowderPatternComponent()):\n", " pdiff = p.GetPowderPatternComponent(i)\n", " if not isinstance(pdiff, PowderPatternDiffraction):\n", " continue\n", " pdiff.SetReflectionProfilePar(ReflectionProfileType.PROFILE_PSEUDO_VOIGT, 0.00001)\n", "\n", "p.UpdateDisplay()\n" ] }, { "cell_type": "markdown", "id": "13", "metadata": {}, "source": [ "### Fit profile, step 2\n", "Refine only constant width, zero, Eta Gaussian/Voigt mix, and a, b, c parameters " ] }, { "cell_type": "code", "execution_count": 8, "id": "14", "metadata": {}, "outputs": [], "source": [ "lsq = LSQ()\n", "lsq.SetRefinedObj(p, 0, True, True)\n", "lsq.PrepareRefParList(True)\n", "# lsq.GetCompiledRefinedObj().Print()\n", "lsqr = lsq.GetCompiledRefinedObj()\n", "\n", "lsqr.FixAllPar()\n", "# lsqr.Print()\n", "# print(lsq.ChiSquare())\n", "lsq.SetParIsFixed(refpartype_scattdata_scale, False)\n", "for par in [\"W\", \"Zero\", \"Eta0\", \"a\", \"b\", \"c\"]:\n", " for i in range(p.GetNbPowderPatternComponent()):\n", " # This is a KLUDGE - we need this because parameter names are\n", " # unique, and thus \"U\" gets renamed to \"U~\", \"U~~\" in case of 2,3 phases.. \n", " lsq.SetParIsFixed(par + \"~\"*i, False)\n", "lsq.SafeRefine(nbCycle=10, useLevenbergMarquardt=True, silent=True)\n", "\n", "p.UpdateDisplay()\n" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "### Fit profile, final\n", "Refine more parameters, and fit the scale factor" ] }, { "cell_type": "code", "execution_count": 9, "id": "16", "metadata": {}, "outputs": [], "source": [ "lsqr.FixAllPar()\n", "# lsqr.Print()\n", "# print(lsq.ChiSquare())\n", "lsq.SetParIsFixed(refpartype_scattdata_scale, False)\n", "for par in [\"U\", \"V\", \"W\", \"Zero\", \"Eta0\", \"Eta1\", \"a\", \"b\", \"c\", \"2ThetaDispl\", \"2ThetaTransp\"]:\n", " for i in range(p.GetNbPowderPatternComponent()):\n", " # This is a KLUDGE - we need this because parameter names are\n", " # unique, and thus \"U\" gets renamed to \"U~\", \"U~~\" in case of 2,3 phases.. \n", " lsq.SetParIsFixed(par + \"~\"*i, False)\n", "lsq.SafeRefine(nbCycle=10, useLevenbergMarquardt=True, silent=True)\n", "\n", "p.FitScaleFactorForRw()\n", "p.UpdateDisplay()\n" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "### Compute weight percentages\n", "This uses the formula: \n", "$w_i = \\frac{S_iZ_iM_iV_i}{\\Sigma_iS_iZ_iM_iV_i}$\n", "\n", "where:\n", "* $w_i$ is the weight fraction of crystalline phase i\n", "* $S_i$ its scale factor in the Rietveld refinement\n", "* $Z_i$ the multiplicity of the formula in the unit cell\n", "* $M_i$ the crystal formula's molecular weight\n", "* $V_i$ the unit cell volume for the phase\n", "\n", "This assumes that the structure is known (and thus that the\n", "CIF files are correct), and that we know all present phases.\n", "\n", "The obtained numbers can be compared to:\n", "https://www.iucr.org/__data/iucr/powder/QARR/results.htm" ] }, { "cell_type": "code", "execution_count": 10, "id": "18", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Weight percentages:\n", "Aluminium oxide - $-alpha: 2.42%\n", " Zincite: 4.18%\n", " Calcium fluoride: 93.40%\n" ] } ], "source": [ "# TODO: check if the method is correctly applied, notably\n", "# is there a shortcut in ObjCryst++ so that the calculated\n", "# structure factor sometimes skips a factor 2 (centrosymmetry\n", "# or other centering factors which allow to avoid a direct sum)\n", "\n", "w = p.qpa(verbose=True)\n" ] } ], "metadata": { "kernelspec": { "display_name": "objcryst", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.5" } }, "nbformat": 4, "nbformat_minor": 5 }